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ABSTRACT 



Context. The strength and direction of migration of embedded low mass planets depends on the disc's thermodynamic state. It has 
been shown that, in discs where the viscous heating is balanced by radiative transport, the migration can be directed outwards, a 
process which extends the lifetime of growing planetary embryos. 

Aims. In this paper we investigate the influence of opacity and stellar irradiation on the disc thermodynamics. We focus on equilibrium 
discs, which have no net mass flux. Utilizing the resulting disc structure, we determine the regions of outward migration in the disc. 
Methods. We perform two-dimensional numerical simulations of equilibrium discs with viscous heating, radiative cooling and stellar 
irradiation. We use the explicit/implicit hydrodynamical code NIRVANA that includes a full tensor viscosity and stellar irradiation, 
as well as a two temperature solver that includes radiation transport in the flux-limited diffusion approximation. The migration of 
embedded planets is studied by using torque formulae. 

Results. In the constant opacity case, our code reproduces the analytical results corresponding to a black-body disc: the stellar 



irradiation dominates in the outer regions - leading to flaring (H/r 



,.2/7 



) - while the viscous heating dominates close to the star. In 



particular, we find that the inner edge of the disc should not be significantly puffed-up by the stellar irradiation. If the opacity depends 
on the local density and temperature, the structure of the disc is different, and several bumps in the aspect ratio Hj r appear, due to 
transitions between different opacity regimes. The bumps in the disc structure are very important, as they can shield the outer disc 
from stellar irradiation. 

Conclusions. Stellar irradiation is an important factor for determining the disc structure and has dramatic consequences for the 
migration of embedded planets. Compared to discs with only viscous heating and radiative cooling, a stellar irradiated disc features a 
much smaller region of outward migration for a range of planetary masses. This suggests that the region where the formation of giant 
planet cores takes place is smaller, which in turn might lead to a shorter growth phase. 

Key words, accretion discs - planet formation - hydrodynamics - radiative transport - planet disc interactions - stellar irradiation 



1. Introduction 

In accretion discs around young stars planetary embryos are born 
and grow through collisions to form bigger and bigger objects. 
The largest protoplanets constitute the cores of giant planets. 
Embedded planets interact with the gas in the disc and can move 
through th e disc. This process depends on the disc's physical 
properties dWardlll997h . in particular the local density and tem- 
perature and the gradients thereof. 

In a real disc, not only viscous heating and radiative cool- 
ing determine the disc's thermal structure, but also stellar irra- 
diation. In the outer par ts of the disc the heating f rom the star 
can keep the disc flared ( IChiang & Goldreichll 19971) . in contrast 
to the disc profile without stellar irradiation. The effect of stel- 
lar irradiation on the disc structure has largel y been inve s tigate d 
utilizing a 1D+1D num erical approach (e.g. iBell et all (1 19971) : 
iDullemond et al.l (1200 ll) ) with the goal of fitting the spectral en- 
ergy distributions (SEDs) of observed discs 

The migration of a low-mass planet embedded in a 
fully radiative gaseous disc can be significantly different 
from migration in an isothermal or purely adiabatic disc 
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2010). While all authors agree that radiation transport can slow 
the rate of inward migration, there is still a lack of consensus 
whether the direction of migration can be o utward and, if so, in 
which part of the disc dBitsch & Klevll201 ll) . Part of this confu- 
sion may be due to the sensitivity of th e direction and magnitude 
of migration on local disc parameters dPaardekooper et al.ll20ldl 
iMasset & Casoh1l2010t IPaardekooper et al.ll201 ll). including, for 
exam ple, the radial disc temperature gradient dAvliffe & Bate! 
1201 ll) . 

In all previous works on planetary migration, the disc struc- 
ture is determined by viscous heating and cooling. The equilib- 
rium disc structure is th us determined by the disc mass and the 
magnitude of viscosity dBitsch & KlevluOlll) . Including stellar 
irradiation can lead to a flared disc pro file in the outer parts 
of the disc dDullemond & Dominikll2004l) . The resulting change 
in the disc structure in the outer parts can have a huge effect 
on the migration of embedded planets. As the aspect ratio of 
the disc H/r changes, so does the temperature profile, which 
in turn determines the strength of the entropy gradient. The en- 
tropy gradient determines the magnitude of the corotation torque 
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(iBaruteau & Massell2008l) . which constitutes a significant frac- 
tion of the total torque, and hence determines the strength and 
direction of migration. 

In this paper we investigate the structure of the disc due to 
changes in opacity and the effects of stellar irradiation using a 2D 
(r - z) disc model. In Section[2]we describe the energy equation 
including stellar heating and the opacities used. We then com- 
pare test studies with analytical calculations of the disc structure 
in the case of constant opacity in Section [3] followed by studies 
with non constant opacity in Section|4] The differences between 
discs including stellar heating and those with only viscous heat- 
ing are described in Section l4~3l We then apply torque formulae 
to estimate the torque acting on embedded planets for these disc 
structures in Section [5] In Section [6] the summary and conclu- 
sions are presented. 

2. Methods 

The protoplanetary disc is treated in this study as a two- 
dimensional (2D) non-self-gravitating gas whose motion is de- 
scribed by the Navier-Stokes equations. In this paper we focus 
on the axisymmetric structure of the disc and compute a verti- 
cal slice (in the r — z plane) utilizing 2D spherical coordinates 
(r-ff). Turbulence in discs is thought to be driven by magneto- 
hydrodynamical instabilities (Bal bus & HawlevHl998l) but here 
we treat viscosity utilizing either a constant kinematic viscosity 
or an a-viscosity dShakura & Sunvaevll973l) . The dissipative ef- 
fects can then be described via the standard visc ous stress-tensor 
approach (eg. lMihalas & Weibel Mihalasll984l) . We also include 
the irradiation from the central star, which will be described in 
detail in section. [27X1 For that purpose we modified and substan- 
tially extended an existing multi-dimensional hydrody namical 
code Nirvana dZiegler & Yorkdll997tlKlev et alj|200lb . 

The radiative energy associated with viscous heating and 
stellar irradiation is then diffused through the disc and emitted 
from its surfaces. To describe this pr ocess we utilize the flux- 
limite d diffusion approximation (FLD. lLevermore & Pomraningl 
1 198 lb . an approximation that allows us to transition from the op- 
tically thick mid-plane to the thin regions near the disc's surface. 

The hydrodynamical equatio ns solved in the c ode have al- 
ready been described in detail dKlev et al.l 120091) . so we fo- 
cus here on our newly implemented two-temperature approach 
in Section 12.11 and how to solve the coupled equations in 
AppendixIBI 

2.1. Energy equation 

In the two-temperature approach (featuring the radiative energy 
density Er and the thermal energy density e) the evo lution equa- 
tions for the thermal and radiation energy read (lKlevlll989t 
iDobbs-Dixon et al.ll201(A ICommercon et al.(l201 lb : 



+ V-F = P k p (T,P)[B(T)-cE r ] 



(1) 



8Er 
dt 

^ + (u-V)e = -PV-u- P k p (T,P)[B(T)-cE r ] + Q + +S. 

Here, Er is independently evolved from the thermal component 
with B(T) = 4oT 4 (<x being the StefanBoltzmann constant and 
T the temperature of the gas). F denotes the radiative flux. From 
eq.Q] one can calculate the coupling time-scale T coup = l/(p*r) 
with the density p, the opacity k, and c the speed of light. The 
coupling time-scale defines the coupling between B(T)/c and 
Er. If pK is high, the coupling time is small, indicating that Er re- 
laxes quickly towards B(T)/c, meaning Er sb 4crT 4 /c (optically 



thick limit). If pK is low, the coupling time is large, and Er and 
B(T)/c can decouple (optically thin limit). B{T)-cEr represents 
the exchange of energy between the thermal and radiative com- 
ponents through emission and absorption of low energy photons. 
Kp is the Planck opacity (see section l2~3l ). The thermal energy e 
is given by e — pc v T (c v being the specific heat at constant vol- 
ume which we hold fixed and uniform), and u = (u r , ug, u v ) the 
gas velocity. P is the gas pressure, Q + the viscous dissipation 
function and S the contribution from the stellar heating. 

The radiative flux, using flux-limited diffusion (FLD), can be 
written as 



Ac 

F = VE R , 

PKr 



(2) 



where kr is the Rosseland mean opacity and A the flux lim- 
iter. Please note that Kp and kr are different opacities. More de- 
tails concerning the opacities used here are given in Section [2~3l 
Using FLD allows us to solve for stable accretion disc models 
that cover several vertical pre ssure scale heights. We use here 
the FLD approach des cribed inlLevermore & Pomraningl d!981l) 
with the flux-limiter of Kiev ( 1989) given by 



A = 



3+ V9+10R 2 
10 



1CW+9+ V81 + 180S 



for R < 2.0 
for R > 2.0 



where 



R = 



1 \VE R \ 



(3) 



(4) 



PKr Er 

These are the same specifications as in iKlev et al.1 (120091) . who 
used a one-temperature approach with E r = a r T 4 , where a r - 
4-cr/c is the radiation constant. 

For our purposes the physical extent of the star is small 
with respect to the disc extension so we approximate it as a 
point source. In spherical coordinates, the stellar radiation is thus 
propagated along the radial direction only. With this method, the 
code simultaneously treats the stellar heating component via ray- 
tracing and subsequent re-radiation via FLD. The stellar heating 
density S (i.e. the energy deposited per unit time and volume), 
received by a grid cell i of width Ar is given by: 



1 



■ e 



-p«*Ar 



Ar 



with F+ 



Ria-Tl 



(5) 



where P* is the stellar radius and the stellar surface temper- 
ature. The factor e~ T < expresses how much stellar irradiation has 
been absorbed in the grid before it hits the radial grid cell ;, with 
Tj being the optical depth integrated radially towards grid cell i. 
The factor 1 - e~ pK * Ar describes how much stellar irradiation is 
absorbed in the actual grid cell, with /e* bei ng the optical opac- 
ity, c alculated using the stellar spectrum dDobbs-Dixon et all 
120101) . For more information regarding the opacities, please see 
Section. l2~3l In the optically thin limit (pK+Ar < 1.0), we make 
the following approximation: 



S = F+e T, pK+ 



(6) 



More details concerning the stellar heating function are given in 
Appendix lAl 

The coupled energy equations (eq. [TJ have to be solved 
simultaneously and we follow here the approach outlined in 
ICommercon et al.l (1201 ll) . The advection (u ■ V)e and the com- 
pressional heating (-PV-u) terms are solved in separate routines. 
Appendix IB1 describes our approach to simultaneously solving 
the coupled energy equations (eq.QTJi. 
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2.2. Physical setup 



2.2.2. Stellar heating 



2.2.1. Computational Parameters 

For computational reasons, the inner disc (0.0416AU < r < 
1.04AU) and outer disc (1.04AU < r < 50.96AU) are treated 
in two different sets of simulations. Both sets of simulations are 
2D discs in the r-9 direction (radial and vertical) with 384 x 32 
(inner disc) or 384x64 (outer disc) active grid cells. The opening 
angle used by the grid, i.e. the range of 9, differs for the inner and 
outer disc simulations. In the inner disc 83° < 9 < 90°, while in 
the outer disc 70° < 9 < 90°, where 9 is the colatitude measured 
from the z-axis and 9 = 90° is the mid-plane. 

The reason for introducing two sets of simulations for the in- 
ner and outer disc is numerical. The time step for the inner disc, 
due to the small inner radius, is very small, significantly increas- 
ing the computation time. It is therefore not feasible to simulate 
the whole disc in one simulation. Additionally, the opening an- 
gle of the numerical disc needs to increase to larger values at 
larger radial distances as the initial disc structure is flared (H/ r 
increasing with f). Setting an opening angle too large (e.g. 20°) 
for the inner disc results in a collapse of the time step due to very 
low densities in the upper regions of the disc. Test runs using a 
density floor showed the same behaviour. Therefore, the two sets 
of simulations utilize different opening angles. The simulations 
then can be attached to each other at the respective boundaries of 
the numerical grids, resulting in a continuous profile. The trans- 
fer of the stellar heating is described more precisely in section 
|2"T21 

We assume the upper and lower parts of the disc are iden- 
tical and therefore apply symmetric boundary conditions at the 
disc's mid-plane, i.e. at 9„ mx = 90°. At the top of the disc (at 
Omin) we set the radiation energy to Er = orT 4 with T = 3.0K 
(the temperature of the interstellar medium). In this way the disc 
will always be cooled by the upper boundary and all the heating 
created in the disc (viscous and stellar) can be transported out- 
wards, as the boundary radiation energy is always lower than in 
mid-pla ne. This type of c ooling at the top of the disc has been 
used in lKlev et al.N2009l) . a nd is also adopted in SPH simula- 
tions of fully radiative discs (lAvliffe & Batell20Tll) . In the radial 
direction we use reflecti ng boundary cond itions for all disc pa- 
rameters as described in lKlev et al.l d2009l) . 

We start the simulations with an initial surface density 
profile of S = E WlAUr 1/2 , where £ = 1000g/cm 2 or 
So - 3000g/cm 2 , wh ich is approximately the value used in 
IWeidenschillingl (I1977I) at 1AU. The initial aspect ratio H/r oc 
r 2 ^ 1 indicates a flar ed disc, following the "flaring disc principle" 
dDullemondll2002l) : if the disc can flare, it will. Self-shadowed 
discs are discs that cannot be made to flare, which may occur for 
discs that are initialized with constant H/r. 

We use either a constant kinematic viscosity of v = 
10 15 c m 2 /s or an a prescription wi th different values for a, fol- 
lowing [Shakura & Sunvaevl (1 19731) with v = ac 2 /Q.K, where c s 
is the mid-plane sound speed and Q.k the Keplerian orbital fre- 
quency. In case of a-viscosity we adopt a vertically constant vis- 
cosity utilizing the mid-plane values. MHD simulations indicate 
that the turbulent stresses do not strongly change with height 
dFlaig et al.ll2012l) . For discs with constant viscosity, the initial 
surface density profile is the equilibrium profile, which cancels 
viscous transport. For ff-discs the profile will evolve until a new 
equilibrium profile is achieved under the influence of reflecting 
boundary conditions. In this way, the mass inside the disc is con- 
served. This case is what we call equilibrium disc. 



The star at the centre of our grid has = M Q , T+ = 5656K 
and R+ = 3.0R o , corresponding to a typical protostar. In our 
simulations the inner discs starts at r = 0.041 6 AU, close to the 
corotation radius of the star. 

Let us first assume that the inner edge of the disc is sharply 
cut off and it cools through that edge as a blackbody. If all of the 
stellar irradiation is absorbed at this edge, the received energy 
per time is 



o-TiRl 



H 



S = 2nr min 2HF* = Anr min H^-± = AtktTIrI 



(7) 



where 2nr m i„2H is the surface of the inner face of the disc. r m! „ 
denotes the radius of the inner edge of the disc. The blackbody 
cooling of the surface of the inner edge is given by 



Q = -4nr mi „Ho-T 4 D , 



(8) 



with To being the disc's temperature. In equilibrium (heating is 
equal to cooling) we get 



S + Q = 



-4nr min HcrT 4 D = 



T 4 R 2 = T 4 r 2 (9") 
1 * n * 1 D'min- W 



Using vertical hydrostatic equilibrium, 



H Y GM* n 



nun J ' mm 



r min ^ 



(10) 



with fi being the mean molecular weight and ft the gas constant, 
we obtain 



TiRl 



H 

' min 

H 

• mil 



\" G 4 Mjf, 4 2 
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T^RlK 4 2 
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C l JVJ 4 



= 0.0516 £ 
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3R e 



1/4 



) l/2 ( 
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J/4 



I 



5600KJ 



Eq.fTTI indicates that the aspect ratio of the inner rim scales with 
the grid truncation radius (r ml „) as r° ( 25 ; decreasing r mm results 
in a decreased H/r. Therefore, regardless of our choice of r m! „, 
provided that the chosen value of r m( „ is smaller than the radius 
of the innermost bump in H/r due to opacity transitions (see 
Section 0), we are guaranteed that the disc interior to this point 
will not have a larger aspect ratio than that at the chosen trunca- 
tion radius. This has important consequences for the irradiation 
received by this first cell, as it guarantees this first cell will not 
be in the shadow of a higher H/r region interior to it. Though the 
simulation naturally captures self-shadowing in the simulated re- 
gions of the disc, eq.fTTI shows this is not relevant for the region 
interior to r„„„. Therefore, if r m ,„ is larger than the physical trun- 
cation radius of the disc, we can mimic the absorption of stellar 
irradiation due to the inner disc by simply reducing the stellar 
irradiation by a factor of e~ T " m " before it hits the first active grid 
cell. In our case T mner corresponds to the optical depth of the in- 
ner two ghost cells of the numerical grid, that are located inside 
of r m j„. These ghost cells have the same properties (e.g. p, k, Ar) 
as the first active cell. By choosing the values of the ghost cells 
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for e~ Timer , we follow the radial profile of the disc, making Ti„ ner 
consistent with the rest of the simulation. 

As some of the stellar irradiation is absorbed by the inner 
disc, the outer disc will not receive stellar irradiation in the op- 
tically thick mid-plane region. From the simulations of the inner 
disc, we can measure to what angle from mid-plane in the disc 
the stellar irradiation is absorbed, & a bs, by looking at the distri- 
bution of stellar heating in the inner disc. This angle i? fl /„ varies 
with the disc mass and viscosity. In the outer disc simulation, 
the disc will only receive stellar heating for an opening angle §„ 
larger than § a b S (keep in mind that i? = 9 - 90° = represent the 
mid-plane of the disc): 



F* = 



0.0 



for 
for 



> &abs 
< &abs 



(12) 



For the simulations of the outer disc, we reduce the stellar irra- 
diation on the top layers by e~ T " m " as we did for the simulations 
of the inner disc. This has the effect that the transition from zero 
stellar irradiation to stellar irradiation is smoothed out a little 
bit. This, of course, implies that the simulations of the inner disc 
have to be done before the simulations of the outer disc to obtain 
the correct § a b s . 

In principle one could radially integrate the stellar flux from 
the inner boundary of the inner simulation to the outer boundary 
to see what remains. This would then act as the starting stellar 
heating of the simulations of the outer grid. However, as the sim- 
ulations of inner and outer disc have different opening angles, 
this could only be done up to the opening angle of the inner disc. 
Also, as the inner parts are (radially) optically thick, most of the 
stellar irradiation is absorbed in the first cells, resulting in a zero 
stellar irradiation at the outer boundary of the inner disc simu- 
lations. Only the very top layers of the inner disc simulations 
still receive stellar irradiation at their outer boundary, which can 
be mimicked by the introduction of # £ ,/„. In Appendix IA. 31 we 
present test simulations showing that the inner and outer disc 
simulations match perfectly with our simple prescription of § a b s . 

2.3. Opacities 

In the previous sections we have introduced three different opac- 
ities, the Planck mean opacity Kp, the optical opacity and the 
Rosseland mea n opacity kr. In principle the opacities can be 
calculated as in iDobbs-Dixon et al.l (l2010t) . The Rosseland and 
Planck mean opacities are defined in the usual manner, so the 
local Planck mean opacity for the low-energy photon group is 
given by 



Kp(T,p) = 



J K v , ns (T,p)B v {T)dv 



(13) 



j B Y (T)dv 

while the optical opacity is given by the high-energy photon 
group 

j /c, vlJ (r,p)7 v (r*)c/v 



(14) 



The frequency dependent opacity is given by k v and the subscript 
ns indicates scattering processes, which are neglected when cal- 
culating wavelength dependent opacities. In principle, the spec- 
trum of impinging radiation can differ from a blackbody, but for 
our purposes we set J V (T+) = B V (T+). Finally, the Rosseland 
mean for the low-energy group is given by 



Kr{T, P ) = F g STr) 



0B r (T) 
dT 



dv 



f 



dT 



dv 



(15) 



p = 10 *g/cm^ 
p = 10" 10 g/cm 3 
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Fig. 1. Opacity from iBeil & Or] (11994 on a logarithmic scale. 
Different colours indicate different densities. The bumps in the 
profile originate from transitions in the opacity regime. This plot 
features kr and Kp as we assume in this paper that Kp = Kp. 



The wavelength-de pendent opacities include t he effect of scat- 
tering (subscript s). IDobbs-Dixon et al 1 (l2010h state that the ra- 
tio between and Kp is about a factor of 10. They calculated 
the opacities directly from wavelength-dependent opacities us- 
ing gaseous, solar-composition opacities. The net result of their 
work is that the photospheres for the stellar irradiation and the 
cooling are physically disconnected. The upper layers of the disc 
(where stellar irradiation plays the largest role) is largely de- 
pleted of dust, particularly at a later stage of evolution, when 
planets are fo rming. Thus th e opac ity due to grains is reduced 
(as is done in lAvliffe & Batel (120091) ). In our first set of simula- 
tions, we leave the Rosseland and Planck opacities constant at 
1 .0cm 2 /g, but the optical opacity is reduced relative to the other 
used opacities by a factor of 10, so that = 0.1cm 2 /g. 

For si mulations with var ying opacities, we follow the opac- 
ity law of IBell & Linl (Il994l) . The opacity depends primarily on 
temperature, but also slightly on density, as can be seen in Fig.Q] 
The plot shows the Rosseland mean and Planck mean opacity. If 
divided by a factor of 10, the opacities in Fig. Q] reflect the op- 
tical opacity. The opacity profile shows several bumps, which 
are caused by transitions in the opacity regime. For example, the 
decrease in opacity around 100K is associated with the melting 
of ice grains, that dominate the opacity at cooler temperatures. 
This opacity law is applied for Kp and kr, although we point out 
here the Rosseland mean opacity and the Planck opacity are not 
the same in reality. However, for the first approximations they 
are quite similar. Again, the optical opacity k+ is 0. 1 of the other 
opacities. 

A change in the opacity, will ultimately lead to a change 
in the discs hydrodynamical structure. Higher values of the 
Rosseland mean opa city lead to a ho t ter and vertically more ex- 
tended disc structure. ISemenov et al.l (120031) also provided opac- 
ity tables for the Rosseland mean opacity and the Planck mean 
opacity, however the differences between the two do not seem to 
be very significant, so we use here the same value for Kp and kr. 
The influence of different opacity laws on the disc structure will 
be investigated in more detail in a future paper, where we will 
also include a more realistic prescription for the optical opacity. 
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Fig. 2. Aspect ratio H/r of the stellar irradiated inner disc cal- 
culated using the mid-plane values of c s for discs with constant 
opacity, k = lcm 2 /g and constant kinematic viscosity. Starting 
from a flared disc profile the discs are evolved until they reach 
thermal equilibrium. 

3. Constant opacity discs 

In this section we focus on the structure of discs with a constant 
opacity, Kp — kr — 1.0cm 2 /g and = 0.1cm 2 /g. We split the 
simulations into an inner disc and an outer disc so we do not 
have the same time step limitation due to the inner edge of the 
computational grid. We present in this section viscous simula- 
tions with both constant viscosity of v = 10 15 cm 2 /s and with 
a viscosity with a = 0.001,0.004,0.008, where v = ac 2 /Q, as 
well as simulations of non-viscous discs. All simulations include 
stellar heating and the discs all feature So = 1000g/ cm 2 . S(r) is 
constant for the simulations with constant viscosity, while S(r) 
changes shape in the discs with a viscosity until an equilibrium 
is reached. 

3.1. Inner disc 

3.1.1. constant-viscosity discs 

We expect the structure of the inner disc to be dominated by 
viscosity. The viscous heating is stronger in the inner parts of 
the disc compared to the stellar heating because Q + oc r~ 3 and 
S oc r~ 2 (see Equation [16]). In Fig. [2] we present the aspect ratio 
H/r of discs with constant viscosity. H is defined as H — c s /CIk, 
where c s is the isothermal sound speed c s - JPjp computed in 
the mid-plane of the disc. The discs are evolved from the initially 
flared structure until they reach thermal equilibrium. 

The non-viscous disc (blue line in Fig. |2j features a puffed- 
up inner rim due to stellar irradiation, while the disc behind 
is shadowed by the rim leading to a vertical contraction and 
a drop of H/r in the outer parts of the simulations. In the- 
ory, if the puffed up inner rim is high enoug h, the whole outer 
disc can be shielded from stellar irradiation dDullemondll2002t 
iDullemond & Dominild 120041) . However, we find that beyond 
0.15AU the discs aspect ratio increases with r, even if the disc 
remains in the shadow of the inner rim. This is because the heat 
propagates outwards from the inner rim. Therefore the tempera- 
ture drops with distance as T — T$r~P. As the disc is in hydro- 
static equilibrium, the aspect ratio increases if /3 < 1 .0. 

For the viscous disc (red line in Fig. |2]i only a very small 
puffed up inner rim is visible, with a marginally increased H/r 



near the inner edge of the disc. The otherwise constant aspect 
ratio is typical of a viscously passive disc, which we show below 
analytically. 

For a purely passive disc (without stellar irradiation), con- 
sider the viscous heating Q + and the radiative cooling Q (black- 
body cooling) for an annulus with size A = InrSr 

q 

Q = AnrSraT* /T efl and Q + = -ZvQ 2 2nrdr . (16) 

8 

Please be aware that the disc can cool from both sides, hence the 
factor 4 in the cooling term. The effective optical depth is given 
by T e ff = 0.5at£. By using the H/r profile for discs in hydrostatic 
equilibrium 

_(H\ 2 GM^ 

Td -\7) 

and with r e ft = 0.52/c we can equate heating and cooling to find 
the aspect ratio 

H l a \l/8 

r ~ y32 ji? G^M\a J r (U > 
H i v \l/4/ \(l-2s)/8 

7 = 0-051 (Toofc) (raj) 

where £ = I.o(r/lA\jy s . Here we have assumed both con- 
stant opacity (a: = lcm 2 /g) and constant kinematic viscosity. 
Therefore, for discs with S = So(r/lAU) -0 5 the result is a disc 
with a constant aspect ratio, which matches quite well with our 
numerical simulations. 

In total the aspect ratio of the viscous disc is higher com- 
pared to the non-viscous disc. This implies that the innermost 
bump of the disc in the non-viscous disc (which is due to stellar 
irradiation) does not influence the structure of the viscous disc, 
as the viscous heating dominates. 

3.1.2. a-viscosity discs 

For a-viscosity discs we also expect the viscosity to dominate 
over stellar irradiation in the inner disc. The corresponding H/r 
profiles are shown in Fig. [3] The black (dotted) lines in the plot 
show the analytical expectation that neglects stellar irradiation. 
In fact the aspect ratio of the different a discs vary with the value 
of viscosity. As the viscosity changes with r, the gradient of the 
surface density changes as well, influencing the aspect ratio. 

For a purely passive disc with a-viscosity, the aspect ratio 
can be computed in the same way as for a constant viscosity 
disc. We equate viscous heating Q + and blackbody cooling Q : 

q 

Q = 4nrSro-T 4 /T eB and Q + = -SvQ 2 2nr6r , (18) 

8 

where v = ac s H = ac 2 /Q.x = aH 2 y£l K . With the relation of the 
hydrostatic equilibrium, the heating of the disc is given by 

Q + = -I.a'R-TQ, K 2nr6r . (19) 

8 yU 

By equating heating and cooling, the aspect ratio of the disc can 
be determined: 

— _Mv2 7f« 1 A 1/6 1/4-W3 f2m 

— = 0.0572 y^Qi) y iooo g °/cm 2 j (tau; ' 
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Fig. 3. Aspect ratio H/r of the stellar irradiated inner disc calcu- 
lated using the mid-plane values for discs with <r-viscosity and 
constant opacities. The black lines indicate the analytical expec- 
tations for the aspect ratio (see text). 



which clearly indicates that for a large enough gradient of the 
surface density s, the aspect ratio can decrease with r. This is 
exactly what happens for our a = 0.008 simulation for which we 
measure s = 0.85. As the inner disc is dominated by viscosity, 
this calculation gives an estimate for the aspect ratio of the discs 
that is in good agreement with our simulation (black lines in 
Fig.0. 

3.2. Outer disc 



Theoretical calculations of IChiang & Goldreichl (fl997h have 
shown that the outer disc is flared with H/r oc r 2 ^ . This can 
easily be shown by equating the stellar heating with the cool- 
ing at the top of the disc. In our simulations for constant opacity 
we find exactly this behaviour, as can be seen in Fig. [4] where 
we display the aspect ratio of discs with viscosity and without. 
For the simulations with viscous heating an absorbing angle of 
$abs = 6.5° has been used, while for the non-viscous simula- 
tions d-gbs = 5.0° has been used, both taken from the results of 
the inner disc simulations. The aspect ratio of the outer viscous 
disc matches perfectly with the inner disc (Fig. [2J, while there 
is some small difference for the non-viscous disc. In the non- 
viscous case the shadowed region behind the innermost rim (in 
the inner disc simulation) seems not to be captured perfectly by 
the continuing outer disc simulation. The non-viscous disc sim- 
ulations show a slight mismatch in the Hj r profile between the 
inner disc and the outer disc. This mismatch cannot be seen in 
the viscous disc simulations, where H/r matches very well. In 
reality, accretion discs are viscous, so we are confident that our 
code reproduces the transition between inner and outer disc very 
well in that case, which we will focus on in section[4] 

It can be seen clearly that both discs show a flared profile in 
the outer regions of the disc, which follows the predicted r 2 ^ 1 
profile. In the outer regions, both disc profiles are nearly the 
same as stellar heating dominates over viscous heating. In the 
inner regions of the simulation, small differences between the 
viscous and non-viscous discs arise, which are due to the impor- 
tance of viscous heating. 

In most 2D simulations of locally isothermal discs (in r - <f> 
plane), a constant Hj r profile is assumed, which corresponds to 
a constant opacity disc without stellar irradiation. However, in 




Fig. 4. Aspect ratio H/r of the stellar irradiated outer disc cal- 
culated using mid-plane values for discs with constant opacity 
k - lcm 2 /g and constant viscosity. Starting from a flared disc 
profile the discs are evolved until they reach thermal equilib- 
rium. These simulations differ from the ones presented in Fig. [2] 
only by the opening angle and radial extent of the disc. 



radiative discs with temperature dependent opacities the aspect 
ratio will depend on radius. This is discussed below. 



4. Disc structure with temperature dependent 
opacity 

In realistic accretion discs, the opacity is not constant. The opac- 
ity depends on the temperature and density (Fig.[T]i. As the tem- 
perature inside an accretion disc drops from a few 1000K in the 
inner parts to a few 10K in the outer parts, we expect the opac- 
ity to vary by orders of magnitude as it was already stated in 
ISemenov et al.l d2003l) and as can also be seen in Fig. Q] As the 
inner parts of the disc are dominated by viscous heating and as 
accretion discs are viscously driven, we focus here only on sim- 
ulations including viscous heating, featuring both constant vis- 
cosity and a viscosity. All simulations shown here also include 
stellar irradiation. 



4. 1. Inner disc structure 
4.1.1. Constant viscosity 

In this subsection, all simulated discs have a constant viscosity of 
v = 10 15 cm 2 /s. In Fig. [5]the aspect ratio and the corresponding 
opacity kr of the inner disc with So = 1000g/cm 2 is displayed. 
Compared to the initial profile the aspect ratio increases, due to 
viscous heating. This effect is also be seen for a constant opacity 
(Fig. |2]i. The drop behind the puffed up inner edge was not visi- 
ble in the simulations with constant opacity (red line in Fig. [2j. 
As can be seen in Fig.Q] the large features in the opacity profile 
at high temperature are important for the inner disc. 

At the inner edge of the grid (r - 0.041 6 AU) the tempera- 
ture drops from ss 10000K to about 3000K in a few grid cells, 
which is associated with a large drop in opacity (top in Fig. |5j. 
As the opacity drops, the cooling rate increases, as it is oc 1/kr. 
An increased cooling also means that the temperature drops even 
further. As the temperature is linked to H/r, a drop in tempera- 
ture translates into a drop in aspect ratio. The drop of H/r behind 
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Fig. 5. Aspect ratio H/r of the stellar irradiated inner disc (bot- 
tom) calculated using the mid-plane values for So = 1000g/cm 2 
discs with non constant opacity as in Fig.[T]and the correspond- 
ing opacity values for kr (top). Starting from a flared disc profile 
the discs are evolved until they reach thermal equilibrium. 



the puffed up inner edge is directly related to a drop in opacity, 
as we show in Appendix lA.2l 

If the inner edge were be shifted towards smaller radii, the 
height of the puffed up inner region would not increase indef- 
initely, as the opacity does not continue to increase, but rather 
levels off at high temperatures (see Fig. [T}. 

As soon as the temperature drops to lower values, the opacity 
increases again (Fig.Q]and top of Fig. [5]), which in turn increases 
H/r. At around r = 0.83AU we observe another bump in the H/r 
profile. Here the temperature crosses the « 1000K bump in the 
opacity profile, which leads to a decrease of H/r for r > 0.83AU. 

We can clearly relate the transition regions of the opacity to 
transition regions of the H/r profile of the disc. These changes 
in the disc's profile are of crucial importance for irradiated discs. 
The bump at r — 0.83AU is higher than the bump at the in- 
ner edge, indicating that stellar irradiation will be absorbed at a 
larger height from mid-plane, which means in turn that the disc 
will only be heated above a larger angle &ai, s . As the inner re- 
gion of the disc is dominated by viscous heating, the height and 
location of the bumps in the disc will change with disc mass and 
the amount of viscous heating. In contrast, stellar irradiation is 
relatively unimportant in the inner disc and does not change sig- 
nificantly with disc mass. 




0.5 0.6 
r[AU] 

Fig. 6. Aspect ratio H/ r of the inner disc (with viscous and stellar 
irradiation) calculated using the mid-plane values for discs with 
a- viscosity and non constant opacity as in Fig.Q] Starting from a 
flared disc profile the discs are evolved until they reach thermal 
equilibrium. 



4.1.2. a-viscosity 

The amount of viscous heating has a crucial influence on the disc 
structure in the interior of the disc. We now focus on discs with 
different a viscosities, with v = ac 2 JQ.K, where we use c s in the 
mid-plane of the disc. The results of simulations of the inner disc 
with different a and Zo values are shown in Fig. [6] 

As expected the aspect ratio of the disc increases with in- 
creasing a parameter due to the increased viscous heating. Also, 
the aspect ratio increases as we increase So, as more mass inside 
the disc results in a larger heating of the disc. This has the effect 
that the peaks in the H/ r profile are shifted towards larger radii, 
as the temperature increases and therefore the transitions in the 
opacity are reached at larger distances from the star. 

The simulation of constant viscosity (with Eo = 1000g/cm 2 
in Fig. [5]l most closely resembles the a = 0.004 simulation. 
However, the value that matches best might be around a as 
0.0055. 

With increasing disc mass and viscosity an increase of the 
height of the peaks in the discs is visible. These bumps in the 
H/r profile are not directly related to the absorption angle § a b s 
of the disc, but it is safe to say that a higher bump in the disc 
results in a larger absorption angle with important effects on the 
structure of the outer disc, as discussed below. 

4.2. Outer disc 

4.2.1 . Constant viscosity 

We have shown in Section [3] that for discs with constant opac- 
ity the discs are fla r ed wi th H/r oc r 2 ^ 7 in agreement with 
IChiang & Goldreichl dl997h . However, for non constant opaci- 
ties we expect another feature in the disc's structure as the tem- 
perature crosses 100K (Fig.[T). The aspect ratio of the outer disc 
structure and the corresponding opacity Kr for two different sur- 
face density values are presented in Fig. [7] 

At r = 1 .04 AU the aspect ratio of the S = 1000g/cm 2 outer 
disc does not exactly match the aspect ratio computed for the 
inner disc (red line in Fig.|5}. In this region the disc is shadowed 
from direct stellar irradiation and any radial transfer of energy 
relies on the re-radiated radiative flux. The outer simulation does 
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Fig. 7. Aspect ratio H/r of the stellar irradiated outer disc (bot- 
tom) calculated using the mid-plane values for discs with con- 
stant viscosity and non constant opacity as in Fig.Q]and the cor- 
responding opacity values for kr (top). Starting from a flared 
disc profile the discs are evolved until they reach thermal equi- 
librium. The different colours indicate different So-values. The 
red line represents the outer disc simulation corresponding to 
the inner disc simulation shown in Fig. [5] 

not include a source of the re-radiated flux at the inner boundary, 
leading to the small difference in the inner disc. 

The disc with So = 1000g/cm 2 features another bump in the 
H/r profile at r ~ 7.3AU. At this point in the disc the tempera- 
ture crosses T ~ 100K, which is the transition region for melted 
ice grains to solid ice grains. The increase of H/r for r < 7.3AU 
is related to the increase of opacity for 200K > T > 100K, which 
can clearly be seen in the top of Fig. [7] The minima and max- 
ima of opacity and H/r are not exactly at the same location, 
but are slightly shifted. However, a clear trend is still visible. 
Interestingly, for r > 7.3AU the H/r profile monotonically de- 
creases and the profile follows that of an non-irradiated disc due 
to the self-shadowing of the disc. 

The fact that the disc is not flared in the outer parts of the disc 
is a combination of the value of the surface density and of the 
opacity. In the constant opacity scenario, the outer parts of the 
disc were flared (Fig. [4j> for the So = 1000g/cm 2 disc. However, 
in reality for lower temperatures (T < 100K), the opacity drops 
to very small values (see top in Fig. [7}. Please keep in mind 
that k+ = 0. Ikr. As the opacity drops, so does the optical depth 
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Fig. 8. Aspect ratio H/r of the outer disc (with or-viscosity and 
stellar irradiation) calculated using the mid-plane values for 
discs with non constant opacity as in Fig. Q] Starting from a 
flared disc profile the discs are evolved until they reach thermal 
equilibrium. These simulations differ from the ones presented in 
Fig. |6]only by the opening angle and radial extent of the disc. 



t = px+Ar of each grid cell, which means that less and less 
stellar irradiation is absorbed in the outer regions of the disc, as 
the absorption is proportional to t (eq.|6]l. 

In order to keep the disc flared, the upper layers in the outer 
parts of the disc have to absorb stellar irradiation. As the opacity 
drops in this region of the disc, an increase in density would keep 
the disc flared, as t = /e+pAr. Indeed, this can be seen in Fig. [7] 
where the blue line indicates a disc with So = 3000g/cm 2 . The 
outer parts are flared and follow the 2 /7th profile. Also interest- 
ing to note is that the bumps in the inner region of the disc in 
the Hj r profile are shifted to larger distances from the star com- 
pared to the lower mass disc (red in Fig. |7J. This is due to the 
fact that the higher mass disc produces more viscous heating and 
is therefore warmer. We compare the disc with So = 3000g/cm 2 
in great detail with a disc that is only heated through viscosity in 
Sectionl4~3l 

4.2.2. a-viscosity 

The bumps of the inner disc can shield the outer disc from stellar 
irradiation implying the outer disc structure is influenced by the 
inner disc structure, which, in turn, is determined by viscosity 
and disc mass. In Fig.[8]the aspect ratio for discs with different 
ff-viscosity and disc mass are displayed. 

The So = 1000g/cm 2 discs are not flared for any viscosity. 
The outer disc is just too thin to absorb stellar irradiation effec- 
tively and viscosity is negligible there, as was seen the constant 
viscosity disc (Fig. [7J. However, we do see the innermost bump 
in the disc rise as the viscosity increases. 

For the So = 3000g/cm 2 disc with a = 0.001 the equilibrium 
disc structure features a flared disc profile, which is quite similar 
to that of a disc with constant opacity (Fig. [7]). However, for 
a > 0.004, the discs are not flared any more. As the viscosity 
increases the inner bumps of the disc structure and block stellar 
irradiation from going through to the outer parts of the disc. This 
effect increases with increasing viscosity. 

In order to have a disc flared in the outer regions, the disc 
must not only have the right amount of mass, but also the vis- 
cosity of the disc should not be too large, as the viscously domi- 
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Fig. 9. Aspect ratio H/r of the outer disc calculated using the 
mid-plane values for discs with constant viscosity and non con- 
stant opacity as in Fig. [JJ Starting from a flared disc profile the 
discs are evolved until they reach thermal equilibrium. The plot 
features one disc with and one without stellar heating; both discs 
have a constant viscosity. 



nated inner disc can block stellar irradiation from the outer disc. 
Therefore, simulations of the inner disc are always needed to de- 
termine the structure of the outer disc correctly. Whether a disc is 
shadowed or not, is determined through viscosity and disc mass. 



4.3. Passive disc 

We now focus on the difference between stellar irradiated discs 
and discs without stellar irradiation. For this purpose we refer to 
the case of constant viscosity v and So = 3000g/cm 2 , as it was 
already shown in Fig. [7] 

In Fig. [9] the two aspect ratios of the disc with So = 
3000g/cm 2 are displayed. In the inner parts, the aspect ratios for 
both discs are nearly identical, as this region is dominated by 
viscous heating and not by stellar irradiation. Beyond « 20AU 
the stellar irradiated disc is flared, while the disc without stellar 
irradiation collapses to small H/r at large radii. The 100K bump 
at r ~ 9AU is more prominent in the stellar irradiated disc and 
also shifted a little bit towards larger distances compared to the 
non-stellar irradiated disc, as the disc receives more heat which 
shifts the bump. 

The 2D density map of the discs is shown in Fig. [TO] The 
density profile of the disc without stellar heating does not extend 
to the same height, because we only use 83° < 9 < 90° for 
this disc. We reduced the vertical extent of the grid in this case 
because of the smaller disc thickness. In the disc with stellar 
irradiation, on the other hand, the opening angle of the grid is 
larger, as we need the disc to absorb the stellar irradiation in the 
top layers. 

In Fig. [10] the H(r) line is indicated in black. In both cases 
the H(r) line rises, however, only for the stellar irradiated disc 
does H/r increase (Fig.|9]l. On the other hand, the bumps, which 
are clearly visible in the H/r profiles are not visible in the H(r) 
profile, because the increase of H(r) is too strong to notice the 
small bumps. 
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Fig. 10. Density map (in log-scale) of discs with stellar heating 
(top) and without stellar heating (bottom). Both discs undergo 
viscous heating and have So = 3000g/cm 2 . The black dotted line 
indicates H(r) and the red solid line represents the Ti = 1.0 line 
integrated from the top of the disc. 



Of interest is also the red line in Fig. [TO] which refers to the 
t; = 1.0 line integrated from the top (infinity) of the disc. We 
define in this case 



ti(z) = I pK R dz' , 

U DO 



(21) 



where z is the vertical thickness of the disc and the integration 
is performed on lines of constant radius r. This line is roughly 
connected to the flux F (eq. [2]i of the disc and therefore gives 
a rough estimate of the location of the photosphere where the 
disc is cooling. One can clearly see the bumps in the t/ = 1.0 
line for the stellar irradiated disc. These bumps are at the same 
location as the ones in the H/ r profile, which are not visible in 
this figure. As t depends on kr and p, it is no surprise to see a 
drop of density just above the t; = 1.0 line around r « 15.6AU. 
The bumps of the r; line can not be seen in the non-irradiated 
disc. This line just levels off for large distances to the star, where 
the disc must contract from its initial state in order to maintain 
its energy loss to the upper and lower layers. 

One of the most important indicators of the disc structure is 
the temperature of the disc, as the temperature is related to the 
entropy. Due to the stellar irradiation we expect the upper layers 
of the disc to be heated from the star, while the mid-plane re- 
gions far away from the star are heated vertically from the upper 
layers. The temperature in mid-plane might therefore be slightly 
smaller compared to the top layers. The temperature is displayed 
in Fig. Q~TJ Again we over plot the line for t/ = 1.0 (red), H(r) 
(black) and additionally in blue the t* = 1 .0 line, which we de- 
fine as 



T*(r) 



pK+dr , 

Jo 



(22) 
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Fig. 11. Temperature map (log-scale) of a disc with stellar and 
viscous heating. Indicated in red is the 77 = 1.0 line, blue the 
r+ = 1.0 line and in black the H(r) line. The disc has So = 
3000g/cm 2 . 

where the integration is performed along lines of constant 9. 
The t+ line represents where in the disc the stellar heating is 
absorbed effectively. This line does not match the t/ line at 
10AU< r < 25AU, where t* > 17. This means the stellar heating 
of the disc is effective higher above the mid-plane compared to 
the cooling, which results in a higher temperature at the upper 
layers of the disc, as can be seen clearly in Fig.[TT] 

The upper layers of the disc are heated by stellar irradia- 
tion, while the disc's interior is also heated through viscosity. 
However, the heating effect of viscosity is diminished in the up- 
per layers of the disc, as the density is much lower compared to 
the mid-plane of the disc (see Fig.lTOt. In the region of the bump 
in optical depth (10.4AU < r < 23.4AU), the effect of stellar 
heating is reduced as the opacity is reduced, which leads to less 
absorption of stellar irradiation. In this region the upper layers of 
the disc are much hotter compared the regions of the disc right 
below. This temperature inversion in vertical direction of the disc 
could have interesting implications to planets on inclined orbits 
in this part of the disc, as the migration speed is dependent on 
the gradient of temperature. 

In the outer regions of the disc, viscosity is negligible and 
the disc is heated solely by stellar irradiation. The heat is then 
transferred from the upper regions to mid-plane, which causes 
the temperature to show only little vertical dependence. 

Observations indicate that protoplanetary discs have an out- 
wards flared profile, which is consistent with our stellar irradi- 
ated model. The differences in the disc structure between the two 
models are dramatic, and we therefore recommend using stellar 
irradiated disc models for the outer parts of the disc. However, 
as we have seen in Fig. Q and Fig. [8] stellar irradiated discs can 
be self shadowed if the disc mass is small or the viscosity large. 
In both cases the inner bumps of the discs are so large that they 
block stellar irradiation from propagating to the outer parts of the 
disc. For low disc mass models, passive discs follow the same 
H/ r profile as stellar irradiated discs, because the discs are opti- 
cally thin and can not maintain the flared disc profile. 

5. Implications to planet migration 

In order to estimate what is the difference between the torque act- 
ing on planets in discs wit h and without stellar irradia tion, we ap- 
ply the torque formula of IPaardekooper et al.l (1201 ll) to the discs 
described in Section 14.31 These discs feature So = 3000g/cm 2 , 
a constant viscosity v and the opacity law in Fig. [1] The for- 
mula captures the behaviour of the torque caused by Lindblad 



resonances and horseshoe drag on low-mass planets embedded 
in gaseous discs in the presence of viscous and thermal diffu- 
sion. This for mula gives the best m atch to the full 3D radiative 
simulations o flBitsch & Klevid201 ll) . We do not simulate planets 
embedded in stellar irradiated discs in this work, as it is beyond 
the scope of the pa per at this point. 

The formula of IPaardekooper et alJ (1201 ll) is very complex 
and features many details, so that we do not cite every step of the 
formula at this poi nt. The formula is also additionally displayed 
in the Appendix of lBitsch & Klevl d201 ll) . However, we want to 
point out a few key items of the formula. The total torque acting 
on an embedded planet is a composition of its Lindblad torque 
and its corotation torque: 

r M = r L + r c (23) 

The Lindblad torque depends on the gradients of tempera- 
ture T oc and surfa ce density 2 oc r~ ! . It is given in 
IPaardekooper et al.l d201 ll) by 

r r L /r = -2.5- 1.7/3 + OAs and r = (I)' Xprfal , (24) 

where q is the mass ratio between planet and star, h — H/r, l,p 
the surface density of the disc at the planets location and rp the 
distance of the planet to the host star. One can clearly see that a 
change in the gradient of temperature influences the Lindblad 
torque. The same applies to the corotation torque, which is 
strongly dependent on the gradient of entropy, S oc with 
£ = /? - (y - l.0)s. The largest contribution of the corotation 
torque arises from the entropy related horseshoe drag, which is 
given by 

yThs,emlT = 7.9^ . (25) 

y 

The aspect ratio H/ r of the disc with and without stellar ir- 
radiation changes from flared to non flared in the outer parts of 
the disc. This change is related to a change in temperature and 
the local temperature gradient. In Fig. Q~T] the outer parts of the 
disc seem to have a very small radial temperature gradient that 
reduces the effect on the entropy related horseshoe drag. In the 
disc without stellar irradiation, on the other hand, planets can 
migra te outwards to quite la rge distances, depending on the disc 
mass dBitsch & Klevl l20"Tll) . We therefore expect that the zero- 
torque radius for planets will be at a smaller distance to the star 
for stellar irradiated discs compared to discs with only viscous 
heating, as the temperature gradient is larger in the latter case. 

In Fig.[T2]the torque acting on planets with different masses 
in a stellar irradiated disc and a non stellar irradiated disc is dis- 
played. The black lines encircle the regions of outward migra- 
tion. If a planet is inside the black circles it will migrate out- 
wards, but it migrates inwards for the rest of the positions in 
the diagram. At the left side of the black circles, a planet would 
face an unstable torque equilibrium, as the planets would migrate 
away from the line in both directions, while at the right side of 
the circles the planet are in a stable equilibrium as they migrate 
towards the line from both directions. 

We only focus here on planets with a few earth masses, 
as larger planets open gaps in discs and m igrate inwards in 
type-II-migration dLin & Papaloizoulll986allb1;ICrida et alJl2006t 
ICrida & Morbidellil 120071) .' We only display the torque up to 
rp = 18.2AU as the torque acting on planets in the stellar ir- 
radiated disc is negative for larger distances to the star. On the 
other hand, the torque is still positive at that point in the disc 
with only viscous heating. In fact the torque is still positive up 
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Fig. 12. To rque acting on planets with different masses using the 
formula of iPaardekooper et al.l (1201 ll) . The top plot features a 
disc with stellar irradiation and viscous heating (constant v) and 
the bottom plot a disc with just viscous heating. The colour cod- 
ing of the torque has been cut off at T, at = -0.00015, so that 
everything that is black in the figure features a large negative 
torque. The black lines encircle the regions of outward migra- 
tion. 



tors 47 AU. One can also clearly distinguish two different re- 
gions of positive torque, indicating areas of outward migration. 
Finally, the region of outward migration starts at a larger distance 
to the star for the stellar irradiated disc. 

For both disc types the inner region (1.56 < r < 4.7AU) of 
positive torque is nearly identical, due to the fact that the disc 
structure is similar in the inner parts of the disc, as the disc 
there is dominated by viscous heating and not stellar irradia- 
tion. However, as soon as the planetary mass increases above 
AQM Eanh the region of outward migration becomes smaller with 
increasing planetary mass, and the planets migrate inward for all 
distances for Mp > 6QM Earth- 
Planets can migrate outwards as long as the libra- 
tion time-scale is com parable to the cooling time-scale 
(Baruteau & Masset 2008). For large distances to the central star 
dBitsch & KlevfeOl ll) the density is lower in the disc, leading to 
a shorter cooling time-scale that prevents outward migration. 

As the planetary masses increase (Mp > 35M£ flr ,/,) they can 
start to open up partial gaps inside the disc. The positive coro- 
tation torque originates from a regio n very close to the planet 
dKlev & Cr ida 2008;Eix 

et al.ll20 09'). which is then depleted as 
the gap starts to form, thus preventing outward migration. 

The regions of outward migration in the disc relate to the 
///r-profile. As can be seen in Fig. [9] whenever H/r decreases 
in the disc, the planet migrates outwards. When H/r increases 
the planet migrates inwards. This means, whenever we have a 
change in the H/r profile, we change the direction of migration. 
As the bumps and dips in the H/ r profile are influenced by the 
opacity transitions, the opacity determines the migration. As the 
transitions of opacity are shifted away from the star with higher 



viscosity and larger surface density, so is the zero-torque radii in 
the disc. 

To clarify the statement above, keep in mind that the torque 
depends on the surface density, temperature and entropy gradi- 
ent, T = f(s,fi,%). 



Hccr- h 



» T(r) oc r - l - 2b 



GM± H r -2b-l 
r R 



(26) 



which indicates if the disc shows a drop off in H/r, the temper- 
ature gradient increases, which therefore increases the entropy 
gradient as f = /? - (y - l.0)s. Therefore the contribution of the 
entropy related corotation torque is stronger if H/r decreases. In 
the flaring part of the disc, where H/r oc r 2 ^ , the entropy gra- 
dient is smaller as b — -2/7, which leads to f3 = 3/7, which 
reduces the entropy gradient resulting in a negative torque. 

On the other hand, it seems that for planet masses with Mp < 
IMEarth the torque is always negative and that planets migrate 
inwards at all radii. The reduction of the positive torque for low 
mass planets is due to the fact that the horseshoe region narrows, 
thus the horseshoe drag becoming less pronounced, resulting in 
a smaller net-torque acting on the planet. 

For these small planetary masses the inward migration speed 
is about the same for the passive and stellar irradiated disc. It 
also seems that the minimum mass a planet needs to sustain out- 
ward migration is increased in the case of a stellar irradiated 
disc, compared to the non stellar irradiated case. The reason 
for this is a slight change in the entropy gradient of the disc, 
which is steeper in the non stellar irradiated disc. However, we 
believe that even smaller mass planets can migrate outwards as 
the torque formula shows several issues that should be kept in 
mind: 



- The formula bv IPaardekooper et al.l (1201 ll) is dependent on 
the smoothing of the planetary potential, which influences 
the torque. A change in the smooth ing changes the torqu e 
acting on the planet (see Appendix in lBitsch & Klevl d201 ll) ). 
For Fig.[T2]we used a smoothing of 0.4. 

- The formula shows so me differences to full 3D simulations 
dBitsch & Klevl 120 111) , as it was derived from 2D simula- 
tions. 

- 3D radiation-hydrodynamical simulations have shown that 
small mass plane ts with 5ME ar th can still migrate outwards 
dKlev et al.ll2009t) . 

In order to verify these assumptions, high resolution 3D simu- 
lations with embedded small mass planets must be done, as the 
minimal mass for outward migration has not been investigated 
with enough detail. We will address this issue in a further paper. 

However, the torque formula gives an estimate of the 
strength of migration in the disc. The difference in the range of 
outward migration for the two disc types does not change the 
implication of the existence of the zero-torque radius. At zero- 
torque distance, planetary embryos and protoplanets can gather 
and collide and then form larger objects. This is of crucial im- 
portance for the growth of massive planets, where a core of a 
few earth ma sses needs time to a ccrete gas in order to form a gas 
giant planet dPollack et al.ll 19961) . 



6. Summary 

We have investigated the influence of opacity and stellar irradia- 
tion on the structure of protoplanetary accretion discs. Utilizing 
our calculated disc structures we have estimated the torque on 
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embedded planets in discs with and witho ut stellar irradiation 
by usin g the theoretical torque formula from lPaardekooper et al.l 

ooTih . 

Before investigating the influence of opacity on the disc 
structure, we compared non-viscous and viscous discs in the in- 
ner and outer regions of the disc for a constant opacity. In the 
non-viscous disc, a puffed-up inner rim shields the first part of 
the disc from stellar irradiation, while the outer part of the disc 
is flared with H/r oc r 2 ^ 7 , as was sh own in theoretical calcula- 
tions bv lChiang & Goldreichl (Il997l) . In the viscous case, in the 
inner part of the disc, viscous heating dominates over stellar ir- 
radiation. By equating viscous heating with radiative cooling, 
we have shown that the aspect ratio H/r is constant until stellar 
heating dominates. This was con firmed by our simu lations. 

We follow the opacity law of iBell & Lir] ( 1 19941) for our non 
constant opacity discs. As the opacity changes due to transi- 
tions in the opacity law, so does the disc structure. We have seen 
that the bumps and dips in the opacity law reflect the bums and 
dips in the discs structure for stellar irradiated and non irradiated 
discs. If the disc is more massive, the disc produces more viscous 
heating, which results in a higher temperature at larger distances. 
But as the temperature determines the opacity, the bumps in the 
disc are moved outwards with increasing temperature. 

The outer parts of the disc (r > 7.8AU) can be flared if the 
disc mass is high enough. The effect of the disc mass was not 
visible in the constant opacity simulations, as a constant opacity 
leads to a larger optical depth in the upper, less dense regions of 
the disc. The optical depth is crucial in determining how much 
stellar irradiation is absorbed. If the disc is optically thin, the 
stellar irradiation just passes through the disc without any heat- 
ing effects. Therefore the disc absorbs less heat and can not sus- 
tain the flared profile. A minimum mass inside the disc is needed 
in order to keep the discs flared. 

For increasing viscosity, the aspect ratios of the discs in- 
crease. This also means that the bumps due to the opacity in the 
disc become higher. A higher bump in the disc leads to more ab- 
sorbed stellar irradiation in the region of the bump. If the bump 
is high enough, it prevents stellar irradiation from reaching the 
outer disc, resulting in a non-flared disc. Very high viscosity 
therefore leads to self-shadowed discs (Fig. [8]). 

In the inner parts of the disc, where viscous heating domi- 
nates, the disc structures for discs with and without stellar irra- 
diation are similar. In the outer regions, the stellar irradiated disc 
can be flared, which is not the case for the non-stellar irradiated 
disc. In the outer parts of the stellar irradiated disc, where the 
viscous heating is negligible, the disc shows a small dependence 
in T{z), as the heat is transported from the upper layers towards 
mid-plane. 

To estimate the torque acting on planets in discs with stellar 
irradia tion, we apply the torque formula by iPaardekooper et ail 
d20Tlh . where we expect the different disc structure in the outer 
parts to result in different migration scenarios. The discs feature 
two regions of outward migration. One in the inner part of the 
disc (1.56 < r < 4.7 AU) that is nearly identical for the discs with 
and without stellar irradiation, as viscous heating dominates in 
that region of the disc. The second region of outward migration is 
completely different for the two discs. It is much smaller for the 
case of stellar irradiation, as the very outer parts of the disc have 
nearly a constant temperature, which favours inward migration. 
In the case of only viscous heating, outward migration continues 
to much larger distances from the star (r * 47 AU). 

However, important here is that the disc features regions of 
outward migration for different planetary masses. At the corners 
of these regions (zero-torque radius), planetary embryos can ac- 



cumulate and grow until the cores are large enough to accrete 
gas and form gas giant planets. We note that the results from 
this torque formula are only an estimate and real 3D simulations 
with embedded planets need to be done, especially for low mass 
planets, as they all would migrate towards the star according to 
the torque formula. 

As the differences between stellar irradiated discs and pas- 
sive discs are dramatic in the outer regions of the disc, we rec- 
ommend using stellar irradiated disc models for high mass discs. 
Low mass discs are self-shadowed, so that the structure is well 
captured by passive disc models as well. 
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Appendix A: Stellar heating 

A. 1. Calculation of stellar heating in the code 

To compute the amount of stellar heating that is deposited in a 
grid cell, we compute the difference between what arrives and 
what leaves a grid cell: 



f>- r ' - F+e-™ = F*e- T > 1 



= F*e" Ti (l - £>- (T <+'- T <>) = F+e"" (l - e'P""^) , 



(A.l) 
(A.2) 



where i + 1 marks the i + 1th grid cell. The total flux emitted by 
the star is given by L* = 47i7? 2 t cr7^:. The stellar flux per surface 
area on a sphere with radius r is thus 



(A.3) 



With the front area of a grid cell given by 

A = r 2 A^)(cos 6\ - cos 2 ) (A.4) 

we can compute the flux received by a single grid cell 

F+ A = RlcrTlAipicosdi - cos 6 2 ) = F s + , (A.5) 

which then leads to the stellar heating s of a grid cell 

s = F s ^(l-e-"« Ar ) (A.6) 

= R*(rT* Atfcos 6»i - cos 6 2 )e- Ti (l - e - piKihr ) . (A.7) 

As the energy equation (eq. [TJ is written for energy densi- 
ties, we have to divide s by the volume of the grid cell V = 
r 2 ArA^(cos 6>i - cos 82) to get the same units and ultimately to 
get S = s/V as it is used in the code: 



1 _ e -p**Ar 

A~r 



(A.8) 



This result yields a dependence on resolution. However, this 
is only important for the optically thick regions (pK+Ar > 1); 
for optically thin regions we make the approximation (1 - 
e~ pK * &r ) / (Ar) = pK+ . In the inner, optically thick regions of the 
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Fig. A.l. Inner bump of the disc featuring different r ; „,„ for the 
same disc configuration as in Fig. [5] The top panel shows k, 
while H/r is shown in the bottom panel. 

disc the resolution influences the repartition of the stellar heat- 
ing, but it does not influence the global net heating. Resolution 
studies indicate that the height of the puffed up inner edge varies 
by ss 3% when lowering the resolution by 33%. As the disc is 
radially optically thick, the stellar flux is absorbed in the first 
cells; beyond that the simulations match perfectly for all tested 
resolutions. 

A.2. Innermost bumps and relation to opacity 

The innermost bump shown in Fig. [5] and Fig. [6] can be mis- 
taken as a puffed up inner rim caused by stellar irradiation, which 
would be in contradiction to what is stated in eqQT| However, the 
temperature in that region is a few 1000K. At that temperature 
the opacity profile shows a kink (see Fig. [TJ, indicating that the 
disc structure will change. In Fig. [5] the change of H/r can be 
related to a change in opacity. As the opacity drops for slightly 
larger distances than r ml „, so too does H/r. As the opacity deter- 
mines the disc structure, we show in Fig. IA. 1 I simulations of the 
inner disc with different r mm for the same disc configuration as 
in Fig. [5] 

For distances larger than 0.041 6 AU, the opacity drops (top in 
Fig. IA.lt and therefore H/r drops as well (bottom in Fig. IA.U . 
For distances closer to the star, the opacity drops (as already 
indicated by Fig. [TJ because of the higher temperature) and H/r 




r[AU] 

Fig. A.2. Matching of the inner disc simulations with the outer 
disc simulations for an a = 0.008 disc with Eo = 3000g/cm 2 . 

decreases again. This clearly shows that the innermost bump is 
related to a transition in opacity, as the aspect ration of the disc 
drops for smaller and larger distances from 0.0416AU. 

As the opacity for even larger temperatures decreases con- 
tinuously (due to the assumed power-law in the opacity for these 
temperatures) our simulations cover the innermost rim, indicat- 
ing that no additional rim is inside. 

A3. Matching of inner and outer disc 

The models of the inner and outer disc can be attached to each 
other. In order to do so, the regions of the simulations have to 
overlap. In Fig. IA.2I we present the overlapping region of the 
a = 0.008 simulation with "Lq = 3000g/cm 2 , which are actually 
the simulations presented in Fig.|6]and Fig. [8] 

Clearly H/r of the simulations of inner and outer disc match 
very well. We are therefore confident that the approximation of 
6 a bs, as presented in eq.[12] for the stellar heating is sufficient. If 
the matching was done in a region of the disc that is dominated 
by stellar irradiation, the matching might fall into a shadowed 
region of the disc that would then not be captured by the outer 
disc leading to a mismatch of the two simulations. However, if, 
as in our simulations the matching occurs in a viscously domi- 
nated region of the disc, the outer discs inner edge will adapt to 
the outer edge of the inner edge due to the viscous heating. 

Appendix B: Energy equation 

For simplicity the derivation of solving the energy equation 
(eq.[TJ is shown in Cartesian coordinates in the implemented 3D 
algorithm, although we used 2D simulations in the r-6 plane in 
this work. The radiation energy equation is then given by 



E" R +l - El 



Al 



tVe ^R.r+lji ^R.ij.k _ rye ^R.i.j.k n R.r-l,,.k \ 

(+1,;,* Ax i,j.k Ax j 

(Evi+l _pn+\ E"»+l _tvi+l \ 

fVy ^fl.ij+U R.i.j.k _ rVy R.i.jM ^R.i.j-l.k \ 

i.)+U Ay i,j,k Ay J 



1 ,k Ay 

f\Z ZSJJMi R.i,I.k _ pf n Rj.I.k n R,i.j.k-l \ 

i,j,k+l Az i,j,k Az j 

+pK P [4o-(T n+1 ) 4 - CE'^ 1 ] , (B.l) 
where the upper indices (« and n + 1) indicate the time step and 

i>;.u -H'>^ + n - -A ■ 
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D)j,k = h{ D u* + D u-u) . 

are the averaged diffusion coefficients, which are given by 
D u , k = — • (B.2) 

Pi,j,k K R,i,j,k 

The term (J n+1 ) A is non linear and makes the scheme difficult to 
invert, yet it is much easier to solve a linear system implicitly. 
Assu ming that the changes in temperature are small in each time 
step dCommercon et al.ll201 ll) . we can write 



(B.3) 



Now, we need to have T" +l as a function of T", E n R and E" R +l 



in order to solve equation. IB. II We use now the energy density 
equation, where we omit advection and compressional heating, 



'll 

at 



- P k p (T,P)[B{T)-cE r ]+Q + + S 



(B.4) 



With 6 = pc Y T we get 



Al 



— (4cr(r +1 ) 4 - c£" +1 ) + — + . (B.5) 

Cy pcy pcy 



With our approximation for the temperature (eq. IB.3t we find 



r +1 = m + mET 

where 



T" + l2At^(r(T") 4 +™+^- 

Cy v ' pcy pcy 

1 + 16Af^cr(r») 3 

Cy v 7 



and 



m = 



At^c 

cy 



1 + l6At^cr(T") 3 



(B.6) 

(B.7) 
(B.8) 



We can now plug this all into the radiation energy equation and 
solve for E" R +l . We arrive at a matrix equation: 

P\,i,j,kER,i+\,j,k + Pl,i,i,kER,i-l,j,k + P'i,i,i,kER,i,j+\,k + fi4,i,j,kE R jj-\,k 
+fi5,iJ,kE R ,i,j,k+l + P6,iJ,kE R JJ,k-\ +FiJ,kERjj,k = Ri,j,k , 

where the superscript n + 1 has been omitted on the left hand 
side. The matrix elements are given by: 



Written in matrix notation we have 



ME R " +1 = R , 



(B.9) 



which is a matrix that is very similar in nature to the one for 
the one-temperature ene rgy equation and can be solved with the 
same matrix solver as in lKlev et al.l d2009l) . After the matrix iter- 
ation, the new temperature of the system can be calculated with 
eq. IB.6l followed by the computation of e = pc v T. 
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